/*
 * Multimalla.h
 *
 *  Created on: 18/08/2013
 *      Author: jose
 */

#ifndef MULTIMALLA_H_
#define MULTIMALLA_H_
#include <Malla.h>
#include <MetNum.h>
#include <EcDiff.h>
#include <iostream>
#include <fstream>
#include <cmath>

using namespace std;

template <int N>
class Multimalla {
public:
	Multimalla(int i, int j, double x1, double x2, double y1, double y2);
	void Interpolar(int e);
	void Restringir(int e);
	void Corregir(int e);
	void entrada();
	void output_file(int e);
	void Ciclo_V(int N_bajada, int N_mbajo, int N_subida, int N_malto, double Tol);
private:
	Malla m[N];
	MetNum M;
	EcDiff Ec;
};
template <int N>
Multimalla<N>::Multimalla(int i, int j, double x1, double x2, double y1, double y2) {
	float log2i;
	float log2j;

	log2i=ceil(log2f(i-1));
	log2j=ceil(log2f(j-1));
	if(N>=log2i||N>=log2j){
		cout<<"El tamaño de malla elegido no alcanza para el numero de restricciones"<<endl;
		exit(1);
	}
	i=(int)pow(2,log2i)+1;
	j=(int)pow(2,log2j)+1;



		for(int e=(N-1); e>=0;e--) {
			m[e].resize(i,j,x1,x2,y1,y2);
			m[e].SetNivel(e,N);
			i=((i-1)/2)+1;
			j=((j-1)/2)+1;

		}
}

template <int N>
void Multimalla<N>::Interpolar(int e){
	matrix<double> E;
	if ((N-1)>e){
		for (int i=0; i<(m[e].Size_i()-1); i++){
				for(int j=0;j<(m[e].Size_j()-1);j++){
					if(m[e+1].D(2*i,2*j)==0){
						m[e+1].Set_E(2*i,2*j,(m[e].V(i,j)));
						m[e+1].Set_E(2*i,2*j+1,(0.5*(m[e].V(i,j)+m[e].V(i,j+1))));
						m[e+1].Set_E(2*i+1,2*j,(0.5*(m[e].V(i,j)+m[e].V(i+1,j))));
						m[e+1].Set_E(2*i+1,2*j+1,(0.25*(m[e].V(i,j)+m[e].V(i+1,j)+m[e].V(i,j+1)+m[e].V(i+1,j+1))));
					}
				}
			}
	}else{
		cout<<"La funcion interpolar termino con error "<<endl;
		exit(1);
	}
	for (int i=0;i<m[e+1].Size_i();i++){
		for(int j=0;j<m[e+1].Size_j();j++){
			m[e+1].Set_V(i,j,m[e+1].V(i,j)+m[e+1].E(i,j));
		}
	}
}

template <int N>
void Multimalla<N>::Restringir(int e){
	matrix<double> Rtemp(m[e].Size_i(),m[e].Size_j());
	if (e>0){
	if(e==N-1){
		for (int i=0; i<m[e].Size_i(); i++){
			for(int j=0;j<m[e].Size_j();j++){
				if(m[e].D(i,j)==0){
					Rtemp(i,j)=Ec.Get_R(&m[e],i,j);
				}
				else if(m[e].D(i,j)==1)
					Rtemp(i,j)=0;
			}
		}
	} else{
		for (int i=0; i<m[e].Size_i(); i++){
			for(int j=0;j<m[e].Size_j();j++){
				if(m[e].D(i,j)==0){
					Rtemp(i,j)=m[e].R(i,j)-(((m[e].V(i+1,j)-2.0*m[e].V(i,j)+m[e].V(i-1,j))/(pow(m[e].Hx(),2)))+((m[e].V(i,j+1)-2.0*m[e].V(i,j)+m[e].V(i,j-1))/(pow(m[e].Hy(),2))));
				}
				else if(m[e].D(i,j)==1)
					Rtemp(i,j)=0;
			}
		}
	}
	for (int i=0; i<=(m[e-1].Size_i()-1); i++){
			for(int j=0;j<=(m[e-1].Size_j()-1);j++){
				if(m[e].D(2*i,2*j)==0){  //Verificar cuando es necesario tener la matriz D
					if(i==0&&j==0)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i+1,2*j+1)+ 2*(Rtemp(2*i,2*j+1)+ Rtemp(2*i+1,2*j))+ 4*Rtemp(2*i,2*j))));
					else if(i==m[e-1].Size_i()-1&&j==m[e-1].Size_j()-1)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i-1,2*j-1)+ 2*(Rtemp(2*i,2*j-1)+ Rtemp(2*i-1,2*j))+ 4*Rtemp(2*i,2*j))));
					else if (i==0&&j==m[e-1].Size_j()-1)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i+1,2*j-1)+ 2*(Rtemp(2*i,2*j-1)+ Rtemp(2*i+1,2*j))+ 4*Rtemp(2*i,2*j))));
					else if (i==m[e-1].Size_i()-1&&j==0)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i-1,2*j+1)+ 2*(Rtemp(2*i,2*j+1)+ Rtemp(2*i-1,2*j))+ 4*Rtemp(2*i,2*j))));
					else if (i==0)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i+1,2*j-1) + Rtemp(2*i+1,2*j+1)+ 2*(Rtemp(2*i,2*j-1) +Rtemp(2*i,2*j+1)+ Rtemp(2*i+1,2*j))+ 4*Rtemp(2*i,2*j))));
					else if (j==0)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i-1,2*j+1)+ Rtemp(2*i+1,2*j+1)+ 2*(Rtemp(2*i,2*j+1)+ Rtemp(2*i-1,2*j) + Rtemp(2*i+1,2*j))+ 4*Rtemp(2*i,2*j))));
					else if (i==m[e-1].Size_i()-1)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i-1,2*j-1) + Rtemp(2*i-1,2*j+1)+ 2*(Rtemp(2*i,2*j-1) +Rtemp(2*i,2*j+1)+ Rtemp(2*i-1,2*j))+ 4*Rtemp(2*i,2*j))));
					else if (j==m[e-1].Size_j()-1)
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i-1,2*j-1)+ Rtemp(2*i+1,2*j-1)+ 2*(Rtemp(2*i,2*j-1)+ Rtemp(2*i-1,2*j) + Rtemp(2*i+1,2*j))+ 4*Rtemp(2*i,2*j))));
					else {
						m[e-1].Set_R(i,j,(0.0625*(Rtemp(2*i-1,2*j-1) + Rtemp(2*i-1,2*j+1)
							     + Rtemp(2*i+1,2*j-1) + Rtemp(2*i+1,2*j+1)
							     + 2*(Rtemp(2*i,2*j-1) + Rtemp(2*i,2*j+1)
							     + Rtemp(2*i-1,2*j) + Rtemp(2*i+1,2*j))
							     + 4*Rtemp(2*i,2*j))));
					}
				} else {
					m[e-1].Set_R(i,j,0);
					m[e-1].fijar_frontera(i,j);
				}
			}
		}
	}else
		exit(1);

}

template <int N>
void Multimalla<N>::output_file(int e){
	ofstream data ("volt.dat");
		for (int i=0; i<m[e].Size_i(); i++){
			for(int j=0;j<m[e].Size_i();j++){
				data << i << " " << j << " " << m[e].V(i,j) << endl;
			}
				data << endl;
		}
	data.close();
}

template <int N>
void Multimalla<N>::Ciclo_V(int N_bajada, int N_mbajo, int N_subida, int N_malto, double Tol){
	for (int i=N-1;i>0;i--){
//		M.GaussSeidel(&m[i],50,1e-4,Ec,N);
//		Restringir(i);
		M.GaussSeidel(&m[i],N_bajada,1e-4,Ec,N);
		Restringir(i);
	}
	M.GaussSeidel(&m[0],N_mbajo,1e-4,Ec,N);   //esto debería ser solucion directa :s
	for (int i=0;i<(N-1);i++){
		Interpolar(i);
		if (i!=(N-2)){
			M.GaussSeidel(&m[i+1],N_subida,1e-4,Ec,N);
		}
		else{
			M.GaussSeidel(&m[i+1],N_malto,Tol,Ec,N);
		}
	}
}

#endif /* MULTIMALLA_H_ */
